options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # not see parameters str1..., str2,... using regex as exc='[str1,str2,...]' 
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(4,1),mar=c(1,5,1,1))
      drw[,1,] |> plot(type='l',xlab='',ylab=pr)
      drw[,2,] |> plot(type='l',xlab='',ylab=pr)
      drw[,3,] |> plot(type='l',xlab='',ylab=pr)
      drw[,4,] |> plot(type='l',xlab='',ylab=pr)
      par(mar=c(3,5,3,3))
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}

explanatory variable obs.x have noise

xN(x0.sx),yN(b0+b1*x0,s)

ex8-0.stan

normal regression without explanatory variable’s noise

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}

ex9.stan

normal regression with explanatory variable’s noise

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x10;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
  real<lower=0> sx;
  vector[N] x0;
}  
model{
  for(i in 1:N){
    x[i]~normal(x0[i],sx);
    y[i]~normal(b0+b1*x0[i],s);
  }
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x0[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] x1;
  vector[N1] y1;
  for(i in 1:N1){
    x1[i]=normal_rng(x10[i],sx);
    m1[i]=b0+b1*x10[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)

n1=10

#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -88137.4 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
##       33      -37.0164   0.000318117      0.002379           1           1       69    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -37.02
##    b0         7.54
##    b1         1.74
##    s          3.86
##    m0[1]     33.62
##    m0[2]     14.02
##    m0[3]     37.27
##    m0[4]     36.64
##    m0[5]     37.17
##    m0[6]     30.62
##    m0[7]     31.04
##    m0[8]     17.98
##    m0[9]      8.31
##    m0[10]    20.55
##    m0[11]     9.64
##    m0[12]    18.80
##    m0[13]    14.36
##    m0[14]    25.06
##    m0[15]    20.36
##    m0[16]     7.31
##    m0[17]    18.85
##    m0[18]    41.13
##    m0[19]    32.61
##    m0[20]    39.87
##    y0[1]     26.14
##    y0[2]     12.23
##    y0[3]     38.67
##    y0[4]     31.88
##    y0[5]     37.67
##    y0[6]     23.53
##    y0[7]     33.60
##    y0[8]     18.17
##    y0[9]     13.15
##    y0[10]    21.64
##    y0[11]    10.44
##    y0[12]    19.99
##    y0[13]     7.96
##    y0[14]    24.74
##    y0[15]    19.88
##    y0[16]     7.90
##    y0[17]    16.70
##    y0[18]    42.30
##    y0[19]    36.45
##    y0[20]    39.14
##    m1[1]      7.31
##    m1[2]     11.07
##    m1[3]     14.82
##    m1[4]     18.58
##    m1[5]     22.34
##    m1[6]     26.10
##    m1[7]     29.85
##    m1[8]     33.61
##    m1[9]     37.37
##    m1[10]    41.13
##    y1[1]      7.16
##    y1[2]      5.56
##    y1[3]     15.33
##    y1[4]     15.22
##    y1[5]     24.53
##    y1[6]     23.86
##    y1[7]     27.92
##    y1[8]     34.82
##    y1[9]     34.94
##    y1[10]    36.43
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -37.14 -36.86 1.19 1.00 -39.50 -35.84 1.01      625      753
##    b0       7.60   7.55 1.66 1.61   4.83  10.37 1.01      348      421
##    b1       1.74   1.74 0.15 0.15   1.49   1.97 1.01      505      691
##    s        4.35   4.27 0.77 0.76   3.27   5.75 1.00     1498     1310
##    m0[1]   33.62  33.65 1.24 1.18  31.57  35.58 1.00     2579     1537
##    m0[2]   14.07  14.04 1.25 1.21  12.03  16.10 1.01      361      467
##    m0[3]   37.26  37.28 1.46 1.41  34.82  39.56 1.00     1936     1210
##    m0[4]   36.63  36.65 1.42 1.36  34.26  38.87 1.00     2058     1190
##    m0[5]   37.16  37.18 1.46 1.40  34.74  39.45 1.00     1956     1213
##    m0[6]   30.63  30.65 1.09 1.03  28.82  32.34 1.01     2697     1640
##    m0[7]   31.05  31.07 1.11 1.05  29.21  32.79 1.00     2708     1674
##    m0[8]   18.02  18.02 1.05 1.03  16.35  19.75 1.00      418      611
##    m0[9]    8.38   8.33 1.61 1.56   5.69  11.02 1.01      347      420
##    m0[10]  20.59  20.59 0.97 0.94  19.06  22.16 1.00      512      825
##    m0[11]   9.70   9.66 1.52 1.46   7.19  12.19 1.01      348      430
##    m0[12]  18.84  18.84 1.02 1.01  17.23  20.53 1.00      440      693
##    m0[13]  14.41  14.39 1.23 1.20  12.40  16.42 1.01      363      486
##    m0[14]  25.08  25.07 0.94 0.91  23.54  26.60 1.00     1012     1375
##    m0[15]  20.39  20.39 0.98 0.95  18.86  21.97 1.00      502      811
##    m0[16]   7.37   7.33 1.67 1.62   4.59  10.16 1.01      348      421
##    m0[17]  18.88  18.89 1.02 1.01  17.28  20.58 1.00      441      703
##    m0[18]  41.11  41.14 1.72 1.66  38.28  43.86 1.00     1448     1118
##    m0[19]  32.61  32.64 1.19 1.12  30.66  34.48 1.00     2659     1639
##    m0[20]  39.85  39.86 1.64 1.58  37.17  42.44 1.00     1576     1163
##    y0[1]   33.57  33.63 4.58 4.56  25.90  40.92 1.00     1991     2045
##    y0[2]   13.93  14.01 4.51 4.40   6.51  21.32 1.00     1411     1701
##    y0[3]   37.22  37.30 4.60 4.25  29.67  44.45 1.00     2102     2089
##    y0[4]   36.44  36.44 4.66 4.46  28.89  44.42 1.00     1844     1819
##    y0[5]   37.00  37.12 4.73 4.52  29.16  44.30 1.00     2127     1933
##    y0[6]   30.67  30.78 4.64 4.49  23.25  38.15 1.00     2137     2013
##    y0[7]   31.06  30.93 4.54 4.30  23.53  38.60 1.00     2015     1855
##    y0[8]   18.29  18.35 4.65 4.41  10.55  25.80 1.00     1725     1811
##    y0[9]    8.29   8.31 4.85 4.77   0.32  16.00 1.00     1534     1786
##    y0[10]  20.69  20.64 4.50 4.24  13.38  28.20 1.00     1716     1758
##    y0[11]   9.84   9.81 4.64 4.62   2.56  17.41 1.00     1680     1859
##    y0[12]  18.90  18.93 4.48 4.42  11.61  26.23 1.00     1634     1970
##    y0[13]  14.57  14.46 4.58 4.25   7.09  22.20 1.00     1526     1529
##    y0[14]  25.15  25.09 4.59 4.43  17.83  32.57 1.00     1854     1858
##    y0[15]  20.37  20.37 4.44 4.18  13.21  27.69 1.00     1975     2035
##    y0[16]   7.41   7.35 4.74 4.64  -0.40  15.04 1.00     1297     1487
##    y0[17]  18.70  18.92 4.56 4.24  11.14  26.28 1.00     1719     1839
##    y0[18]  41.05  41.00 4.77 4.52  33.11  48.95 1.00     1914     1974
##    y0[19]  32.61  32.60 4.54 4.45  25.50  40.04 1.00     2097     1834
##    y0[20]  39.77  39.90 4.71 4.44  31.53  47.29 1.00     2042     1971
##    m1[1]    7.37   7.33 1.67 1.62   4.59  10.16 1.01      348      421
##    m1[2]   11.12  11.09 1.42 1.38   8.79  13.44 1.01      350      426
##    m1[3]   14.87  14.86 1.20 1.17  12.93  16.85 1.01      367      486
##    m1[4]   18.62  18.62 1.03 1.01  17.00  20.32 1.00      434      693
##    m1[5]   22.37  22.37 0.94 0.93  20.85  23.87 1.00      638     1006
##    m1[6]   26.11  26.11 0.95 0.92  24.55  27.65 1.00     1253     1443
##    m1[7]   29.86  29.87 1.06 1.00  28.11  31.54 1.01     2632     1668
##    m1[8]   33.61  33.64 1.24 1.18  31.56  35.57 1.00     2580     1537
##    m1[9]   37.36  37.38 1.47 1.41  34.92  39.67 1.00     1917     1210
##    m1[10]  41.11  41.14 1.72 1.66  38.28  43.86 1.00     1448     1118
##    y1[1]    7.28   7.21 4.57 4.44  -0.07  14.89 1.00     1302     1651
##    y1[2]   11.17  11.13 4.68 4.61   3.88  18.91 1.00     1304     1698
##    y1[3]   15.06  15.00 4.53 4.15   7.55  22.52 1.00     1499     1554
##    y1[4]   18.71  18.79 4.54 4.18  11.22  26.27 1.00     1770     1912
##    y1[5]   22.41  22.40 4.51 4.43  15.17  29.76 1.00     1756     1834
##    y1[6]   26.17  26.06 4.61 4.52  18.61  33.80 1.00     1895     1788
##    y1[7]   29.66  29.67 4.53 4.25  22.13  37.08 1.00     1959     1847
##    y1[8]   33.69  33.68 4.50 4.10  26.42  41.18 1.00     2094     1727
##    y1[9]   37.40  37.55 4.78 4.32  29.46  45.42 1.00     1881     2000
##    y1[10]  41.19  41.18 4.74 4.51  33.39  48.82 1.00     1882     1892

#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)

mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 1 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 2 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 3 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 3 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 4 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 2 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 4 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 1 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 3 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 1 finished in 0.4 seconds.
## Chain 3 finished in 0.4 seconds.
## Chain 2 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 4 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 2 finished in 0.5 seconds.
## Chain 4 finished in 0.5 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 0.6 seconds.
seeMCMC(mcmc,'[m,x]')
##  variable   mean median    sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -44.46 -45.61 10.56 9.31 -59.89 -25.46 1.05       97      147
##    b0       7.61   7.63  2.02 1.94   4.28  10.73 1.02      444      454
##    b1       1.72   1.73  0.17 0.16   1.44   2.00 1.02      435      349
##    s        2.95   2.98  1.39 1.57   0.74   5.15 1.02      125       85
##    sx       1.77   1.79  0.89 0.88   0.36   3.22 1.02      132      116
##    x0[1]   13.46  13.44  1.46 1.67  11.20  15.66 1.00      429     2495
##    x0[2]    3.78   3.79  1.18 0.98   1.82   5.71 1.01      967     1025
##    x0[3]   16.37  16.48  1.24 1.16  14.30  18.27 1.00      885     1793
##    x0[4]   16.26  16.31  1.17 1.05  14.31  18.09 1.01      996     1930
##    x0[5]   18.74  18.54  1.72 1.81  16.42  21.85 1.02      199      135
##    x0[6]   13.27  13.27  1.14 0.99  11.36  15.02 1.01      747     1579
##    x0[7]   15.69  15.57  1.89 2.25  13.07  18.90 1.02      176      167
##    x0[8]    6.14   6.16  1.20 0.99   4.16   7.99 1.01     1400     1214
##    x0[9]   -1.40  -1.12  1.79 1.85  -4.73   0.91 1.02      214      172
##    x0[10]   6.55   6.64  1.27 1.23   4.41   8.49 1.01      567     1543
##    x0[11]   0.75   0.90  1.25 1.07  -1.43   2.61 1.01      606     1068
##    x0[12]   7.53   7.46  1.35 1.42   5.48   9.65 1.01      275      668
##    x0[13]   5.03   4.98  1.36 1.41   2.97   7.30 1.01      323      446
##    x0[14]   8.73   8.78  1.40 1.48   6.45  10.82 1.01      416     1336
##    x0[15]   9.11   9.05  1.65 1.91   6.81  11.70 1.01      221      763
##    x0[16]  -0.16  -0.10  1.38 1.09  -2.42   1.90 1.02      593      886
##    x0[17]   7.23   7.18  1.23 1.20   5.35   9.21 1.01      480      644
##    x0[18]  20.66  20.44  1.63 1.61  18.42  23.75 1.02      243      139
##    x0[19]  13.54  13.59  1.25 1.22  11.52  15.43 1.00      694     2001
##    x0[20]  17.38  17.43  1.35 1.46  15.21  19.37 1.00      564     1736
##    m0[1]   30.83  30.66  2.62 2.94  27.04  35.20 1.01      299     1955
##    m0[2]   14.19  14.13  2.00 1.68  11.02  17.55 1.01     2631     1856
##    m0[3]   35.80  35.67  2.20 2.13  32.34  39.50 1.01      607     1679
##    m0[4]   35.61  35.42  2.07 1.89  32.31  39.12 1.00      982     2456
##    m0[5]   39.83  40.01  2.70 3.06  35.35  43.65 1.00      348     2469
##    m0[6]   30.47  30.54  1.88 1.61  27.26  33.49 1.01     3531     1951
##    m0[7]   34.59  34.66  2.99 3.70  29.82  38.98 1.01      246     2052
##    m0[8]   18.23  18.23  1.94 1.66  15.09  21.34 1.01     2806     2005
##    m0[9]    5.33   5.16  2.82 3.19   1.37  10.08 1.01      307     1717
##    m0[10]  18.96  18.85  2.18 2.21  15.62  22.66 1.00      586     2779
##    m0[11]   8.99   8.88  2.03 1.83   5.86  12.38 1.01     1677     1477
##    m0[12]  20.61  20.72  2.23 2.34  16.94  24.05 1.01      351     2395
##    m0[13]  16.33  16.52  2.27 2.39  12.54  19.77 1.01      389     1687
##    m0[14]  22.71  22.59  2.38 2.61  19.12  26.52 1.01      354     2123
##    m0[15]  23.31  23.35  2.73 3.16  19.00  27.22 1.01      265     1483
##    m0[16]   7.45   7.50  2.10 1.84   3.96  10.77 1.01     2170     2309
##    m0[17]  20.10  20.16  2.02 1.96  16.81  23.28 1.00      617     2276
##    m0[18]  43.14  43.37  2.50 2.55  38.81  46.74 1.00      566     2375
##    m0[19]  30.95  30.86  2.21 2.23  27.59  34.71 1.01      481     2026
##    m0[20]  37.53  37.29  2.45 2.48  33.94  41.70 1.00      373     1977
##    y0[1]   30.85  30.21  4.12 3.65  25.12  38.32 1.00      759     2570
##    y0[2]   14.18  14.15  3.88 3.17   7.71  20.76 1.00     3469     3390
##    y0[3]   35.81  35.41  4.01 3.29  29.70  42.66 1.01     2016     2266
##    y0[4]   35.53  35.25  3.85 3.08  29.64  42.27 1.00     2742     2643
##    y0[5]   39.83  40.52  4.21 3.68  32.17  45.71 1.00      822     2378
##    y0[6]   30.51  30.52  3.77 3.07  24.15  36.92 1.01     3544     2486
##    y0[7]   34.58  35.19  4.37 4.12  26.62  40.49 1.00      554     2235
##    y0[8]   18.25  18.30  3.74 3.02  12.01  24.48 1.01     3994     2644
##    y0[9]    5.30   4.60  4.32 3.90  -0.60  13.21 1.00      690     2048
##    y0[10]  19.01  18.54  3.95 3.35  13.10  26.06 1.00     1806     2724
##    y0[11]   8.97   8.76  3.81 3.11   2.84  15.49 1.01     3532     3162
##    y0[12]  20.61  21.08  3.93 3.31  13.48  26.38 1.00     1532     3378
##    y0[13]  16.26  16.79  4.07 3.45   8.86  22.17 1.00     1146     2783
##    y0[14]  22.77  22.19  4.08 3.55  16.97  30.25 1.00      974     2469
##    y0[15]  23.36  23.98  4.21 3.70  15.50  29.16 1.00      674     1684
##    y0[16]   7.49   7.49  3.90 3.16   0.99  14.09 1.01     3634     2555
##    y0[17]  20.13  20.40  3.85 3.26  13.66  26.26 1.00     2769     2631
##    y0[18]  43.17  43.75  4.15 3.56  35.68  49.19 1.00     1300     2669
##    y0[19]  30.95  30.46  3.94 3.23  25.29  38.08 1.00     1681     2195
##    y0[20]  37.47  36.89  4.02 3.40  31.56  44.75 1.00     1154     1452
##    m1[1]    7.39   7.40  2.04 1.96   4.03  10.53 1.02      442      454
##    m1[2]   11.11  11.14  1.74 1.66   8.18  13.81 1.02      460      234
##    m1[3]   14.83  14.87  1.46 1.33  12.35  17.17 1.02      490      218
##    m1[4]   18.55  18.58  1.24 1.13  16.45  20.54 1.01      553      213
##    m1[5]   22.27  22.27  1.09 1.01  20.43  24.06 1.01      620      215
##    m1[6]   25.99  25.97  1.07 1.04  24.18  27.76 1.00      630      272
##    m1[7]   29.71  29.69  1.16 1.13  27.84  31.68 1.01      638      611
##    m1[8]   33.43  33.40  1.35 1.29  31.23  35.69 1.01      613     1146
##    m1[9]   37.15  37.16  1.61 1.54  34.38  39.81 1.01      580      261
##    m1[10]  40.87  40.88  1.90 1.81  37.59  44.00 1.01      552      234
##    x1[1]   -0.15  -0.14  1.95 1.43  -3.38   3.01 1.00     3829     3325
##    x1[2]    2.01   2.01  2.03 1.49  -1.20   5.32 1.01     4205     2907
##    x1[3]    4.15   4.15  1.94 1.50   0.95   7.31 1.01     3873     2222
##    x1[4]    6.36   6.34  2.02 1.52   3.14   9.60 1.01     4070     1642
##    x1[5]    8.44   8.45  1.91 1.45   5.35  11.58 1.01     4006     2326
##    x1[6]   10.61  10.62  1.98 1.47   7.34  13.87 1.01     4055     2847
##    x1[7]   12.75  12.80  1.98 1.48   9.44  15.95 1.01     4242     2669
##    x1[8]   14.99  15.01  1.94 1.41  11.85  18.13 1.01     4156     2757
##    x1[9]   17.12  17.14  1.97 1.49  13.91  20.32 1.01     3811     2789
##    x1[10]  19.32  19.28  1.99 1.50  16.16  22.61 1.01     4114     2775
##    y1[1]    7.42   7.43  3.81 3.40   1.35  13.55 1.00     1313     2054
##    y1[2]   11.07  11.14  3.68 3.14   5.15  17.03 1.01     1569     2492
##    y1[3]   14.82  14.80  3.48 2.95   9.09  20.67 1.00     1654     1962
##    y1[4]   18.58  18.54  3.47 2.88  12.96  24.26 1.01     2646     2603
##    y1[5]   22.34  22.23  3.48 2.84  16.79  28.36 1.00     2636     2514
##    y1[6]   25.97  25.95  3.33 2.75  20.48  31.60 1.01     3566     3052
##    y1[7]   29.67  29.56  3.42 2.82  24.18  35.30 1.00     3099     2663
##    y1[8]   33.42  33.35  3.51 2.91  27.81  39.27 1.00     3098     2278
##    y1[9]   37.16  37.10  3.56 3.07  31.42  43.01 1.00     2170     2859
##    y1[10]  40.82  40.75  3.76 3.32  34.81  47.20 1.01     1710     3208

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

outlier

ex10.stan

objective variable have outlier, y~cauchy(b0+b1*x,s)

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~cauchy(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=cauchy_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=cauchy_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)

n1=10

x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -817.234 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       20      -33.1259   0.000203143    0.00056136           1           1       28    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -33.13
##    b0         6.30
##    b1         1.90
##    s          3.18
##    m0[1]     12.01
##    m0[2]     19.76
##    m0[3]     14.92
##    m0[4]     10.58
##    m0[5]      6.38
##    m0[6]      8.18
##    m0[7]     16.78
##    m0[8]     13.25
##    m0[9]     12.50
##    m0[10]     6.91
##    m0[11]     8.62
##    m0[12]    22.06
##    m0[13]    20.27
##    m0[14]    11.03
##    m0[15]    15.72
##    m0[16]     9.03
##    m0[17]    22.80
##    m0[18]    12.08
##    m0[19]    17.88
##    m0[20]    15.65
##    y0[1]     10.05
##    y0[2]     20.11
##    y0[3]     10.98
##    y0[4]      8.94
##    y0[5]      5.57
##    y0[6]      9.06
##    y0[7]     15.88
##    y0[8]     13.73
##    y0[9]      8.03
##    y0[10]     0.33
##    y0[11]     7.52
##    y0[12]    20.30
##    y0[13]    20.03
##    y0[14]    10.39
##    y0[15]    16.68
##    y0[16]    17.01
##    y0[17]    19.68
##    y0[18]    12.49
##    y0[19]    13.15
##    y0[20]    13.93
##    m1[1]      6.38
##    m1[2]      8.21
##    m1[3]     10.03
##    m1[4]     11.86
##    m1[5]     13.68
##    m1[6]     15.51
##    m1[7]     17.33
##    m1[8]     19.15
##    m1[9]     20.98
##    m1[10]    22.80
##    y1[1]      8.43
##    y1[2]      5.43
##    y1[3]      7.78
##    y1[4]      9.01
##    y1[5]     13.08
##    y1[6]      9.33
##    y1[7]     17.45
##    y1[8]     16.43
##    y1[9]     24.82
##    y1[10]    27.25
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -33.55 -33.26 1.27 1.10 -36.15 -32.15 1.00      489      768
##    b0       6.24   6.33 1.52 1.35   3.62   8.76 1.01      280      298
##    b1       1.92   1.91 0.32 0.30   1.40   2.45 1.00      385      486
##    s        3.60   3.51 0.65 0.60   2.71   4.77 1.00     1410     1402
##    m0[1]   11.98  11.98 0.88 0.83  10.50  13.43 1.00      368      569
##    m0[2]   19.79  19.73 1.26 1.19  17.68  21.88 1.00     1292     1424
##    m0[3]   14.91  14.90 0.83 0.78  13.59  16.31 1.00      962      984
##    m0[4]   10.55  10.58 0.99 0.95   8.90  12.19 1.01      308      440
##    m0[5]    6.32   6.42 1.51 1.34   3.73   8.83 1.01      280      303
##    m0[6]    8.13   8.20 1.27 1.17   6.02  10.23 1.01      281      326
##    m0[7]   16.78  16.76 0.93 0.87  15.23  18.33 1.00     1757     1280
##    m0[8]   13.23  13.24 0.82 0.79  11.93  14.61 1.00      498      789
##    m0[9]   12.48  12.49 0.85 0.81  11.06  13.88 1.00      408      654
##    m0[10]   6.85   6.95 1.44 1.29   4.40   9.26 1.01      280      307
##    m0[11]   8.57   8.64 1.21 1.14   6.58  10.59 1.01      283      326
##    m0[12]  22.10  22.04 1.57 1.51  19.48  24.72 1.00      959     1168
##    m0[13]  20.30  20.25 1.32 1.26  18.10  22.52 1.00     1194     1382
##    m0[14]  11.00  11.01 0.95 0.91   9.38  12.59 1.00      322      478
##    m0[15]  15.72  15.70 0.86 0.78  14.33  17.16 1.00     1362     1150
##    m0[16]   8.99   9.05 1.16 1.10   7.07  10.94 1.01      286      353
##    m0[17]  22.85  22.80 1.68 1.63  20.05  25.68 1.00      892     1135
##    m0[18]  12.05  12.06 0.87 0.84  10.57  13.50 1.00      373      585
##    m0[19]  17.90  17.86 1.04 0.96  16.15  19.65 1.00     1714     1165
##    m0[20]  15.65  15.63 0.86 0.78  14.27  17.08 1.00     1330     1150
##    y0[1]   12.02  12.03 3.78 3.61   5.90  18.18 1.00     1736     1998
##    y0[2]   19.66  19.74 3.83 3.68  13.40  25.92 1.00     1673     1654
##    y0[3]   14.95  14.96 3.69 3.48   8.75  20.85 1.00     1772     1853
##    y0[4]   10.40  10.35 3.94 3.85   3.84  16.95 1.00     1733     1691
##    y0[5]    6.38   6.35 4.04 4.07  -0.19  12.98 1.00     1028     1643
##    y0[6]    8.10   8.06 3.80 3.60   1.86  14.48 1.00     1177     1844
##    y0[7]   16.78  16.88 3.81 3.66  10.43  22.67 1.00     1958     1860
##    y0[8]   13.25  13.23 3.80 3.60   6.75  19.33 1.00     1950     1896
##    y0[9]   12.54  12.40 3.72 3.52   6.56  18.93 1.00     1762     1974
##    y0[10]   6.90   6.90 3.86 3.73   0.49  13.18 1.00      957     1642
##    y0[11]   8.45   8.41 3.70 3.52   2.45  14.40 1.00     1531     1807
##    y0[12]  21.96  22.03 3.99 3.90  15.42  28.61 1.00     1826     1740
##    y0[13]  20.44  20.39 3.84 3.80  14.16  26.78 1.00     1874     2060
##    y0[14]  10.83  10.82 3.70 3.66   4.86  16.75 1.00     1342     1884
##    y0[15]  15.74  15.70 3.84 3.77   9.66  21.90 1.00     2070     1974
##    y0[16]   8.91   8.89 3.83 3.67   2.72  15.33 1.00     1387     1815
##    y0[17]  22.80  22.74 4.06 3.89  15.97  29.51 1.00     1834     1942
##    y0[18]  12.14  12.12 3.78 3.62   5.81  18.35 1.00     1339     1695
##    y0[19]  17.99  18.01 3.75 3.36  11.65  24.24 1.00     2035     1800
##    y0[20]  15.53  15.40 3.85 3.74   9.49  21.99 1.00     1998     1920
##    m1[1]    6.32   6.42 1.51 1.34   3.73   8.83 1.01      280      303
##    m1[2]    8.16   8.23 1.26 1.17   6.06  10.26 1.01      281      324
##    m1[3]    9.99  10.03 1.05 0.99   8.28  11.75 1.01      297      423
##    m1[4]   11.83  11.83 0.89 0.85  10.34  13.31 1.00      359      562
##    m1[5]   13.67  13.68 0.81 0.76  12.35  15.03 1.00      576      841
##    m1[6]   15.50  15.48 0.85 0.77  14.13  16.92 1.00     1250     1084
##    m1[7]   17.34  17.31 0.98 0.91  15.70  19.00 1.00     1796     1213
##    m1[8]   19.18  19.13 1.18 1.11  17.20  21.19 1.00     1425     1306
##    m1[9]   21.01  20.95 1.42 1.34  18.66  23.40 1.00     1087     1354
##    m1[10]  22.85  22.80 1.68 1.63  20.05  25.68 1.00      892     1135
##    y1[1]    6.31   6.24 3.92 3.71  -0.14  12.75 1.00     1020     1533
##    y1[2]    8.12   8.25 3.83 3.77   1.74  14.24 1.00     1805     1803
##    y1[3]   10.08  10.06 3.85 3.75   3.70  16.22 1.00     1451     1774
##    y1[4]   11.88  11.88 3.71 3.49   5.78  17.77 1.00     1726     1851
##    y1[5]   13.62  13.63 3.86 3.68   7.16  20.01 1.00     1789     1919
##    y1[6]   15.53  15.59 3.71 3.66   9.49  21.62 1.00     1734     1969
##    y1[7]   17.27  17.29 3.80 3.57  11.07  23.44 1.00     1919     1693
##    y1[8]   19.16  19.23 3.83 3.64  12.91  25.29 1.00     2038     2089
##    y1[9]   20.98  21.06 4.03 3.89  14.29  27.52 1.00     1827     1872
##    y1[10]  22.94  22.92 4.02 4.04  16.48  29.41 1.00     1643     1810

mdl=cmdstan_model('./ex10.stan') 
fn(mdl,data)
## Initial log joint probability = -101.724 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       20       -14.165   0.000738009    0.00153556      0.9772      0.9772       37    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -14.16
##    b0         5.15
##    b1         2.04
##    s          0.49
##    m0[1]     11.28
##    m0[2]     19.59
##    m0[3]     14.40
##    m0[4]      9.75
##    m0[5]      5.25
##    m0[6]      7.17
##    m0[7]     16.39
##    m0[8]     12.61
##    m0[9]     11.81
##    m0[10]     5.80
##    m0[11]     7.64
##    m0[12]    22.05
##    m0[13]    20.14
##    m0[14]    10.23
##    m0[15]    15.25
##    m0[16]     8.09
##    m0[17]    22.85
##    m0[18]    11.35
##    m0[19]    17.58
##    m0[20]    15.19
##    y0[1]     11.57
##    y0[2]     19.59
##    y0[3]     14.87
##    y0[4]      9.36
##    y0[5]     11.07
##    y0[6]      4.19
##    y0[7]     14.81
##    y0[8]     11.15
##    y0[9]     43.30
##    y0[10]     5.98
##    y0[11]     7.72
##    y0[12]    23.50
##    y0[13]    20.67
##    y0[14]    10.31
##    y0[15]    15.04
##    y0[16]     8.02
##    y0[17]    20.20
##    y0[18]    11.54
##    y0[19]    17.53
##    y0[20]    14.52
##    m1[1]      5.25
##    m1[2]      7.20
##    m1[3]      9.16
##    m1[4]     11.11
##    m1[5]     13.07
##    m1[6]     15.03
##    m1[7]     16.98
##    m1[8]     18.94
##    m1[9]     20.90
##    m1[10]    22.85
##    y1[1]     10.20
##    y1[2]      7.50
##    y1[3]      7.32
##    y1[4]     11.18
##    y1[5]     12.88
##    y1[6]     14.78
##    y1[7]     14.36
##    y1[8]     21.08
##    y1[9]     18.71
##    y1[10]    23.10
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median     sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -16.63 -16.32   1.40 1.24 -19.31 -15.04 1.00      438      597
##    b0       5.21   5.19   0.37 0.30   4.65   5.84 1.01      566      482
##    b1       2.03   2.04   0.07 0.06   1.91   2.14 1.01      789      743
##    s        0.74   0.71   0.26 0.24   0.40   1.21 1.01      488      503
##    m0[1]   11.31  11.29   0.24 0.21  10.92  11.71 1.01      711      801
##    m0[2]   19.58  19.58   0.31 0.28  19.07  20.08 1.01     2054     1272
##    m0[3]   14.41  14.41   0.23 0.22  14.04  14.77 1.00     1292     1146
##    m0[4]    9.79   9.77   0.26 0.23   9.37  10.22 1.01      615      648
##    m0[5]    5.30   5.28   0.37 0.30   4.74   5.93 1.01      565      482
##    m0[6]    7.21   7.19   0.32 0.27   6.73   7.76 1.01      562      539
##    m0[7]   16.40  16.40   0.25 0.23  16.00  16.80 1.01     1981     1038
##    m0[8]   12.63  12.63   0.23 0.21  12.26  13.01 1.01      868      919
##    m0[9]   11.83  11.82   0.23 0.21  11.45  12.23 1.01      762      878
##    m0[10]   5.86   5.83   0.35 0.30   5.33   6.46 1.01      562      481
##    m0[11]   7.69   7.67   0.31 0.26   7.22   8.21 1.01      566      544
##    m0[12]  22.03  22.03   0.38 0.33  21.42  22.64 1.00     1746     1245
##    m0[13]  20.13  20.13   0.32 0.29  19.60  20.66 1.01     1982     1312
##    m0[14]  10.26  10.24   0.25 0.23   9.86  10.68 1.01      638      706
##    m0[15]  15.26  15.27   0.23 0.22  14.89  15.64 1.00     1580     1000
##    m0[16]   8.13   8.11   0.30 0.25   7.68   8.63 1.01      572      600
##    m0[17]  22.83  22.83   0.40 0.35  22.18  23.47 1.00     1659     1311
##    m0[18]  11.38  11.37   0.24 0.21  11.00  11.78 1.01      717      801
##    m0[19]  17.58  17.58   0.27 0.24  17.14  18.01 1.01     2146     1075
##    m0[20]  15.20  15.20   0.23 0.22  14.82  15.58 1.00     1558     1004
##    y0[1]   11.88  11.29  22.52 1.09   6.93  16.68 1.00     1826     2015
##    y0[2]   19.65  19.57  12.01 1.08  14.51  24.55 1.00     1921     1430
##    y0[3]   13.51  14.41  36.49 1.06   9.60  18.42 1.00     1829     1988
##    y0[4]   10.31   9.77  31.16 1.14   4.66  14.90 1.00     1888     1878
##    y0[5]    4.62   5.30  20.37 1.12  -0.09  10.11 1.00     1869     1927
##    y0[6]    6.17   7.17  44.66 1.08   2.32  11.65 1.00     2025     2055
##    y0[7]   18.84  16.36  87.97 1.07  11.21  20.85 1.00     1834     1832
##    y0[8]   12.80  12.65  24.53 1.03   7.65  17.02 1.00     1851     2011
##    y0[9]   12.20  11.83  25.40 1.10   6.00  16.77 1.00     2077     1971
##    y0[10]   7.18   5.83  44.20 1.20   0.51  11.61 1.00     1995     2042
##    y0[11]   6.17   7.66  58.92 1.11   2.68  11.82 1.00     1911     1882
##    y0[12]  21.37  22.02  18.48 1.13  17.13  26.02 1.00     1908     1847
##    y0[13]  17.74  20.12 528.57 1.15  15.64  25.02 1.00     1822     1840
##    y0[14]   8.35  10.24  93.64 1.05   6.23  15.08 1.00     1912     1885
##    y0[15]  15.57  15.28  27.93 1.11  10.66  20.60 1.00     1993     2057
##    y0[16]   8.06   8.11  18.39 1.11   3.69  12.47 1.00     2039     1855
##    y0[17]  23.06  22.88  24.58 1.15  18.58  27.78 1.00     1996     2058
##    y0[18]  12.30  11.37  32.48 1.10   6.81  15.95 1.00     1953     1786
##    y0[19]  18.16  17.61  38.75 1.16  13.13  22.46 1.00     2076     2010
##    y0[20]   6.71  15.19 324.96 1.14  10.81  19.74 1.00     2114     1972
##    m1[1]    5.30   5.28   0.37 0.30   4.74   5.93 1.01      565      482
##    m1[2]    7.25   7.23   0.32 0.27   6.77   7.79 1.01      563      558
##    m1[3]    9.20   9.18   0.27 0.24   8.77   9.65 1.01      593      611
##    m1[4]   11.14  11.13   0.24 0.22  10.76  11.55 1.01      699      762
##    m1[5]   13.09  13.09   0.23 0.22  12.73  13.46 1.00      947     1016
##    m1[6]   15.04  15.04   0.23 0.22  14.67  15.42 1.00     1499     1047
##    m1[7]   16.99  16.99   0.26 0.24  16.57  17.41 1.00     2092     1000
##    m1[8]   18.93  18.94   0.30 0.27  18.45  19.42 1.01     2124     1292
##    m1[9]   20.88  20.89   0.34 0.31  20.32  21.44 1.00     1888     1289
##    m1[10]  22.83  22.83   0.40 0.35  22.18  23.47 1.00     1659     1311
##    y1[1]    4.85   5.28  33.30 1.20   0.24  10.33 1.00     1702     1910
##    y1[2]    7.61   7.23  24.64 1.17   2.56  13.01 1.00     1965     1857
##    y1[3]    9.27   9.17  30.60 1.02   5.12  13.33 1.00     1472     1594
##    y1[4]    7.73  11.11 109.98 1.12   6.14  15.54 1.00     1901     1932
##    y1[5]   16.21  13.07 134.29 1.03   8.68  17.60 1.00     1853     2008
##    y1[6]   16.15  15.02  26.11 1.06  10.87  19.56 1.00     1947     1854
##    y1[7]   17.52  17.01  24.32 1.10  12.61  22.50 1.00     1678     1738
##    y1[8]   18.90  18.88  18.82 1.14  13.96  23.59 1.00     2130     1940
##    y1[9]   19.14  20.88 126.44 1.12  16.32  26.22 1.00     2237     2039
##    y1[10]  24.61  22.86  45.36 1.24  18.23  27.63 1.00     1964     1692

censored data

objective variable has NA

(x,y) i=1-N
(x0,y0) i=1-N0
x1 i=1-N1, y1=NA
(x,y)~N((mx,my),(sx2,sy2,sxy))
(x0,y0)~N((mx,my),(sx2,sy2,sxy))
x1~N(mx,sx2)

ex11-0.stan

data{
  int N0;
  array[N0] vector[2] xy;
  int N1;
  vector[N1] x1;
}
parameters{
  vector[2] m;
  cov_matrix[2] s;
}
model{
  target+=multi_normal_lpdf(xy | m, s);
  x1~normal(m[1],s[1,1]^.5);
}
generated quantities{
  vector[2] xy1;
  xy1=multi_normal_rng(m,s);
  real cr;
  cr=s[1,2]/(s[1,1]*s[2,2])^.5;
}
n=30
x=runif(n,0,9)
y=rnorm(n,10+3*x,4)
cor(x,y)
## [1] 0.926
qplot(x,y)

L=4
n0=sum(x>L)
x0=x[x>L]
y0=y[x>L]
x1=x[x<=L]
n1=sum(x<=L)
cor(x0,y0)
## [1] 0.837
qplot(x0,y0)

mdl=cmdstan_model('./ex11-0.stan') 

data=list(N0=n0,xy=matrix(c(x0,y0),ncol=2),N1=n1,x1=x1)

mle=mdl$optimize(data=data)
## Initial log joint probability = -8860.57 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       50      -100.084    0.00102871    0.00225021           1           1       77    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__    -100.08
##    m[1]       4.25
##    m[2]      19.37
##    s[1,1]     6.81
##    s[2,1]    24.16
##    s[1,2]    24.16
##    s[2,2]    96.37
##    xy1[1]     3.39
##    xy1[2]    15.39
##    cr         0.94
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.4 seconds.
## Chain 2 finished in 0.4 seconds.
## Chain 3 finished in 0.4 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 0.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 0.6 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median    sd   mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -96.17 -95.81  1.75  1.63 -99.56 -94.00 1.01      507      815
##    m[1]     4.27   4.28  0.55  0.55   3.37   5.21 1.01      391      388
##    m[2]    19.48  19.62  2.76  2.79  15.04  23.76 1.02      201      115
##    s[1,1]   8.79   8.35  2.69  2.33   5.34  13.89 1.00      814      943
##    s[2,1]  30.80  29.26 11.24  9.99  15.39  51.53 1.02      265      433
##    s[1,2]  30.80  29.26 11.24  9.99  15.39  51.53 1.02      265      433
##    s[2,2] 129.97 118.95 57.31 51.11  56.10 235.03 1.02      210      397
##    xy1[1]   4.28   4.30  2.89  2.81  -0.49   8.89 1.00     1915     1876
##    xy1[2]  19.42  19.91 11.17 10.74   1.22  36.06 1.00     1497      833
##    cr       0.92   0.93  0.06  0.04   0.81   0.97 1.02      476      523

xy=mcmc$draws('xy1')
cor(xy[,,1],xy[,,2])
## [1] 0.889
qplot(xy[,,1],xy[,,2])

objective variable is censored

y i=1-N, y~N(m,s)  
  actual          ya i=1-Na
  lower censored  yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
  upper censored  yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)

cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))

P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)

ex11-1.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  real L;
  int Nl;
  vector[Nl] xl;
  real U;
  int Nu;
  vector[Nu] xu;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
  for(i in 1:Nl)
    target+=normal_lcdf(L | b0+b1*xl[i],s);
  for(i in 1:Nu)
    target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)

xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
mdl=cmdstan_model('./ex8-0.stan')

data=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)

mle=mdl$optimize(data=data)
## Initial log joint probability = -20850.3 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
##       25      -15.3216     0.0011708   0.000321661           1           1       55    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -15.32
##    b0        17.34
##    b1         1.34
##    s          3.33
##    m0[1]     25.91
##    m0[2]     21.57
##    m0[3]     22.06
##    m0[4]     23.04
##    m0[5]     21.61
##    m0[6]     26.18
##    m0[7]     24.32
##    m0[8]     23.91
##    m0[9]     24.67
##    y0[1]     26.97
##    y0[2]     19.61
##    y0[3]     23.36
##    y0[4]     22.76
##    y0[5]     25.05
##    y0[6]     26.71
##    y0[7]     22.03
##    y0[8]     24.08
##    y0[9]     20.02
##    m1[1]     18.89
##    m1[2]     20.06
##    m1[3]     21.23
##    m1[4]     22.39
##    m1[5]     23.56
##    m1[6]     24.73
##    m1[7]     25.89
##    m1[8]     27.06
##    m1[9]     28.22
##    m1[10]    29.39
##    y1[1]     16.50
##    y1[2]     23.37
##    y1[3]     15.09
##    y1[4]     26.74
##    y1[5]     25.46
##    y1[6]     26.43
##    y1[7]     28.54
##    y1[8]     28.90
##    y1[9]     30.68
##    y1[10]    31.85
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m')
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -15.97 -15.55 1.62 1.26 -19.02 -14.31 1.02      356      362
##    b0      18.15  18.17 6.75 6.22   7.99  28.76 1.04      281      213
##    b1       1.19   1.19 1.39 1.22  -1.02   3.25 1.04      310      231
##    s        4.73   4.33 1.82 1.27   2.84   7.83 1.00      843      836
##    m0[1]   25.74  25.77 2.90 2.45  21.05  30.08 1.01      625      521
##    m0[2]   21.90  21.85 2.75 2.39  17.70  26.19 1.02      296      314
##    m0[3]   22.33  22.28 2.37 2.07  18.66  26.07 1.02      321      343
##    m0[4]   23.20  23.14 1.82 1.52  20.33  26.13 1.01      512      707
##    m0[5]   21.93  21.88 2.71 2.35  17.77  26.18 1.02      298      299
##    m0[6]   25.98  26.03 3.14 2.63  20.93  30.66 1.01      571      510
##    m0[7]   24.33  24.35 1.85 1.54  21.33  27.22 1.00     1249      981
##    m0[8]   23.96  23.94 1.73 1.41  21.24  26.72 1.00     1274      972
##    m0[9]   24.64  24.66 2.01 1.66  21.48  27.73 1.00     1100      850
##    y0[1]   25.77  25.87 5.74 5.07  16.74  34.58 1.00     1474     1344
##    y0[2]   21.84  21.65 5.49 4.93  13.50  30.89 1.00      762     1209
##    y0[3]   22.40  22.21 5.80 4.67  13.94  31.49 1.00     1047     1130
##    y0[4]   23.38  23.27 5.51 4.53  15.27  31.62 1.00     1646     1785
##    y0[5]   21.93  21.85 5.80 4.86  13.08  31.22 1.01      879     1301
##    y0[6]   25.76  25.66 6.05 5.04  16.37  35.26 1.00     1268     1215
##    y0[7]   24.30  24.29 5.46 4.76  15.97  33.29 1.00     2003     1541
##    y0[8]   24.14  24.04 5.42 4.68  15.86  32.94 1.00     1778     1535
##    y0[9]   24.61  24.58 5.46 4.82  16.06  33.59 1.00     1859     1891
##    m1[1]   19.52  19.54 5.21 4.73  11.83  27.88 1.04      278      213
##    m1[2]   20.56  20.56 4.09 3.68  14.49  27.09 1.04      278      222
##    m1[3]   21.59  21.52 3.03 2.66  17.00  26.30 1.03      288      260
##    m1[4]   22.62  22.57 2.15 1.89  19.27  25.98 1.01      354      387
##    m1[5]   23.66  23.61 1.71 1.38  20.94  26.42 1.00      851      931
##    m1[6]   24.69  24.71 2.04 1.67  21.44  27.81 1.00     1072      849
##    m1[7]   25.72  25.76 2.89 2.43  21.05  30.04 1.01      628      521
##    m1[8]   26.76  26.83 3.93 3.29  20.42  32.59 1.01      473      453
##    m1[9]   27.79  27.93 5.04 4.26  19.87  35.19 1.01      415      382
##    m1[10]  28.82  28.94 6.19 5.24  19.01  38.01 1.02      385      346
##    y1[1]   19.46  19.53 7.56 6.23   7.83  31.33 1.01      413      602
##    y1[2]   20.51  20.37 6.43 5.52   9.92  30.74 1.01      558      747
##    y1[3]   21.60  21.47 6.02 5.31  12.21  30.95 1.00      853     1221
##    y1[4]   22.69  22.55 5.29 4.67  14.20  31.03 1.00      896     1095
##    y1[5]   23.55  23.57 5.31 4.54  15.22  31.71 1.00     1828     1599
##    y1[6]   24.73  24.62 5.57 4.66  15.94  33.97 1.00     1788     1572
##    y1[7]   25.90  25.98 5.89 4.81  16.41  34.85 1.00     1442     1253
##    y1[8]   26.75  26.74 6.47 5.42  16.71  37.25 1.00      950      927
##    y1[9]   28.02  28.25 7.28 6.18  16.91  39.39 1.01      707      722
##    y1[10]  28.68  28.76 7.93 6.74  15.87  41.41 1.01      541      474

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

data=list(N=n,x=xya$x,y=xya$y,
          L=L,Nl=nl,xl=xyl$x,
          U=U,Nu=nu,xu=xyu$x,
          N1=n1,x1=x1)
mdl=cmdstan_model('./ex11-1.stan') 

mle=mdl$optimize(data=data)
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Initial log joint probability = -333.417 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       25      -21.8357   0.000501952   4.20973e-05      0.9955      0.9955       40    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -21.84
##    b0         5.78
##    b1         3.93
##    s          4.51
##    m0[1]     30.94
##    m0[2]     18.21
##    m0[3]     19.64
##    m0[4]     22.53
##    m0[5]     18.33
##    m0[6]     31.75
##    m0[7]     26.29
##    m0[8]     25.06
##    m0[9]     27.30
##    y0[1]     29.13
##    y0[2]     21.85
##    y0[3]     19.29
##    y0[4]     22.04
##    y0[5]     19.63
##    y0[6]     32.25
##    y0[7]     28.27
##    y0[8]     29.41
##    y0[9]     24.75
##    m1[1]     10.35
##    m1[2]     13.77
##    m1[3]     17.20
##    m1[4]     20.62
##    m1[5]     24.05
##    m1[6]     27.47
##    m1[7]     30.90
##    m1[8]     34.32
##    m1[9]     37.74
##    m1[10]    41.17
##    y1[1]     16.86
##    y1[2]     10.31
##    y1[3]      9.99
##    y1[4]     27.31
##    y1[5]     20.55
##    y1[6]     29.08
##    y1[7]     33.12
##    y1[8]     38.61
##    y1[9]     43.36
##    y1[10]    39.91
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m',ch=T)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -22.08 -21.72 1.48 1.17 -24.85 -20.52 1.02      325      546
##    b0       1.58   2.45 6.69 5.76 -11.81  10.30 1.04      105       94
##    b1       4.88   4.64 1.42 1.20   3.06   7.50 1.04      100      104
##    s        6.42   5.91 2.43 1.81   3.73  10.89 1.02      202      499
##    m0[1]   32.81  32.25 3.44 2.66  28.65  38.72 1.02      204      286
##    m0[2]   17.01  17.41 2.78 2.35  11.60  20.89 1.02      199      271
##    m0[3]   18.79  19.06 2.46 2.10  14.25  22.35 1.02      283      371
##    m0[4]   22.37  22.39 2.09 1.75  19.06  25.47 1.01      701      743
##    m0[5]   17.15  17.53 2.75 2.32  11.83  21.00 1.02      203      293
##    m0[6]   33.82  33.20 3.68 2.83  29.43  40.24 1.02      182      246
##    m0[7]   27.04  26.77 2.32 1.80  24.07  30.69 1.01     1429      716
##    m0[8]   25.52  25.34 2.16 1.74  22.52  28.85 1.01     1410      965
##    m0[9]   28.29  27.98 2.51 1.95  25.18  32.36 1.01      738      619
##    y0[1]   32.73  32.09 7.72 6.28  21.35  45.22 1.01     1541     1105
##    y0[2]   17.28  17.59 7.34 6.48   5.18  28.32 1.01     1174     1388
##    y0[3]   19.01  19.09 7.16 6.17   7.63  30.39 1.00     1685     1257
##    y0[4]   22.56  22.69 7.17 6.32  10.85  34.13 1.00     1714     1165
##    y0[5]   17.42  17.76 7.43 6.35   5.16  28.60 1.01     1102     1042
##    y0[6]   33.86  33.35 7.62 6.52  22.74  46.92 1.01      605      827
##    y0[7]   26.97  26.77 7.16 6.34  15.89  38.60 1.00     1764     1805
##    y0[8]   25.56  25.48 7.18 6.09  14.37  36.87 1.00     1807     1545
##    y0[9]   28.27  28.04 7.12 6.16  17.37  40.04 1.00     1825     1508
##    m1[1]    7.25   8.01 5.14 4.40  -3.19  14.04 1.03      110       88
##    m1[2]   11.50  12.14 4.04 3.49   3.64  16.84 1.03      125      158
##    m1[3]   15.75  16.21 3.04 2.62   9.85  19.84 1.03      171      248
##    m1[4]   20.00  20.18 2.29 1.93  15.94  23.33 1.01      359      433
##    m1[5]   24.25  24.15 2.08 1.70  21.19  27.31 1.01     1186      991
##    m1[6]   28.50  28.17 2.55 1.98  25.36  32.65 1.01      659      564
##    m1[7]   32.75  32.20 3.42 2.64  28.60  38.64 1.02      206      286
##    m1[8]   37.00  36.22 4.47 3.48  31.63  44.99 1.03      146      199
##    m1[9]   41.25  40.36 5.60 4.43  34.57  51.36 1.03      126      124
##    m1[10]  45.50  44.33 6.77 5.43  37.34  57.69 1.03      116      122
##    y1[1]    7.28   8.02 8.66 7.72  -7.80  19.34 1.01      453      388
##    y1[2]   11.53  12.00 7.93 6.73  -2.11  22.87 1.00      803      956
##    y1[3]   15.83  16.08 7.47 6.54   2.89  27.00 1.00      972     1512
##    y1[4]   19.98  20.15 7.53 6.39   7.90  31.54 1.01     1708     1280
##    y1[5]   24.18  24.20 7.33 6.58  12.53  35.73 1.00     1675     1698
##    y1[6]   28.26  27.87 7.52 6.40  17.20  40.63 1.00     1935     1364
##    y1[7]   32.75  32.16 7.67 6.60  21.23  45.45 1.01     1195     1361
##    y1[8]   37.11  36.26 8.27 7.00  25.65  51.35 1.00     1141     1004
##    y1[9]   41.11  40.17 8.99 7.32  28.91  57.03 1.01      387      644
##    y1[10]  45.18  44.09 9.28 7.80  32.69  60.56 1.02      256      437

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

sensitivity/specificity analysis

ex14.stan

estimating sens and spec

data {
  int N;
  array[N] int x;
  array[N] int y;
}
parameters {
  real<lower=0,upper=1> p;
  real<lower=0,upper=1> se;
  real<lower=0,upper=1> sp;
}
model {
  p~uniform(0,1);
  se~uniform(0,1);
  sp~uniform(0,1);
  for (i in 1:N) {
    y[i]~bernoulli(x[i]*se+(1-x[i])*(1-sp));
  }
}
generated quantities {
  real ppv;
  real npv;
  ppv=se*p/((se*p)+((1-p)*(1-sp)));
  npv=(1-p)*sp/(((1-p)*sp)+(p*(1-se)));
}
n=20
data=list(N=n,
          x=sample(0:1,n,replace=T),
          y=sample(0:1,n,replace=T))

mdl=cmdstan_model('./ex14.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -11.4166 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5      -11.2128   0.000275847   1.14736e-05           1           1        8    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -11.21
##      p        0.75
##      se       0.73
##      sp       0.22
##      ppv      0.73
##      npv      0.22
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -17.49 -17.17 1.33 1.11 -20.03 -16.04 1.01      745      996
##      p      0.51   0.52 0.29 0.37   0.07   0.95 1.01      662      632
##      se     0.69   0.71 0.13 0.13   0.47   0.88 1.00     2735     1625
##      sp     0.27   0.25 0.12 0.12   0.09   0.50 1.01     1882     1217
##      ppv    0.50   0.51 0.29 0.38   0.06   0.95 1.01      675      653
##      npv    0.46   0.45 0.30 0.39   0.03   0.95 1.01      748      786

2x2 cross table

Effect occur y=1 with probabilty p01, p11 from each cause x{0,1}
event frequncy nxy of effect y{0,1} by cause x{0,1}
n01~B(n0.,p0)
n11~B(n1.,p1)

n01=n0p0, n00=n0(1-p0)
n11=n1p1, n10=n1(1-p1)

p00=n00/n=n0(1-p0)/n, p01=n01/n=n0p0/n
p10=n10/n=n1(1-p1)/n, p11=n11/n=n1p1/n

Cramer'V  (chi2/n/(min(row,column)-1))^.5
  in 2x2  
  crv =(n11*n00-n10*n01)/(n0.*n1.*n.0*n.1)^.5
      =(n0(1-p0)n1p1-n0p0n1(1-p1))/(n0n1(n0(1-p0)+n1(1-p1))(n0p0n1p1))^.5

kappa coefficient   k=(po-pe)/(1-pe)
  po: Observed agreement (proportion of times both raters agreed)
  pe: Expected agreement under independence
      po=p00+p11
        =(n0(1-p0)+n1p1)/n
      pe=(p00+p01)(p00+p10)(p10+p11)(p01+p11)
        =n0/n*(n0(1-p0)+n1(1-p1))/n*(n0p0+n1p1)/n*n1/n

ex16-1.stan

data {
  int n;
  int n0;
  int n01;
  int n1;
  int n11;
}
parameters {
  real<lower=0,upper=1> p0;
  real<lower=0,upper=1> p1;
}
model {
  n01~binomial(n0,p0);
  n11~binomial(n1,p1);
}
generated quantities {
  real RR;
  RR=p1/p0;
  real OR;
  OR=(p1/(1-p1))/(p0/(1-p0));
}

ex16-2.stan

data {
  int n;
  int n0;
  int n01;
  int n1;
  int n11;
}
parameters {
  real<lower=0,upper=1> p0;
  real<lower=0,upper=1> p1;
}
model {
  n01~binomial(n0,p0);
  n11~binomial(n1,p1);
}
generated quantities {
  real RR;
  RR=p1/p0;
  real OR;
  OR=(p1/(1-p1))/(p0/(1-p0));
  real crv;
  crv=(n0*(1-p0)*n1*p1-n0*p0*n1*(1-p1))/(n0*n1*(n0*(1-p0)+n1*(1-p1))*(n0*p0+n1*p1))^.5;
  real k;
  real po;
  po=(n0*(1-p0)+n1*p1)/n;
  real pe;
  pe=n0/n*(n0*(1-p0)+n1*(1-p1))/n*(n0*p0+n1*p1)/n*n1/n;
  k=(po-pe)/(1-pe);
}
n0=30
n01=rbinom(1,n0,0.3)
n1=30
n11=rbinom(1,n1,0.6)
data=list(n=n0+n1,n0=n0,n01=n01,n1=n1,n11=n11)

mdl=cmdstan_model('./ex16-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -44.3649 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5      -38.8529   0.000121331   1.41282e-05           1           1        8    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -38.85
##      p0       0.30
##      p1       0.57
##      RR       1.89
##      OR       3.05
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -42.82 -42.51 0.97 0.74 -44.73 -41.86 1.01      981     1123
##      p0     0.31   0.31 0.08 0.08   0.18   0.45 1.01     2044     1378
##      p1     0.56   0.56 0.08 0.09   0.42   0.70 1.00     2172     1722
##      RR     1.95   1.83 0.67 0.54   1.11   3.24 1.00     1972     1454
##      OR     3.40   2.93 1.94 1.48   1.21   7.12 1.00     2062     1553

data=list(n=n0+n1,n0=n0,n01=n01,n1=n1,n11=n11)

mdl=cmdstan_model('./ex16-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -47.0484 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5      -38.8529    0.00190612    0.00050074           1           1        8    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -38.85
##      p0       0.30
##      p1       0.57
##      RR       1.89
##      OR       3.05
##      crv      0.27
##      k        0.61
##      po       0.63
##      pe       0.06
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -42.82 -42.51 0.97 0.74 -44.73 -41.86 1.01      981     1123
##      p0     0.31   0.31 0.08 0.08   0.18   0.45 1.01     2044     1378
##      p1     0.56   0.56 0.08 0.09   0.42   0.70 1.00     2172     1722
##      RR     1.95   1.83 0.67 0.54   1.11   3.24 1.00     1972     1454
##      OR     3.40   2.93 1.94 1.48   1.21   7.12 1.00     2062     1553
##      crv    0.26   0.26 0.12 0.12   0.05   0.45 1.00     2082     1633
##      k      0.60   0.60 0.06 0.06   0.49   0.70 1.00     2093     1659
##      po     0.63   0.63 0.06 0.06   0.52   0.72 1.00     2092     1659
##      pe     0.06   0.06 0.00 0.00   0.06   0.06 1.00     1809     1386

point of subjective equality PSE

PSE: 50% threshold for sensing the difference between two stimuli is equal
JND: Just noticeable difference, difference between 50% threshold and 75%

r~B(n,p) #reaction for stimuli
p=1/(1+exp(-(a+b*x)))
x=x1-x0, x0,x1 is strength of stimuli

PSE=-a/b
JND=(log(0.75/0.25)-a)/b-PSE

ex6-3-0.stan

mulit logistic regression

data{
  int N;
  int m;
  vector[N] x;
  array[N] int y;
}
parameters{
  real b0;
  real b1;
}
model{
  y~binomial_logit(m,b0+b1*x);
}
n=20
m=10
x=runif(n,-2,2)
y=rbinom(n,m,1/(1+exp(-(-2+3*x))))

glm(y/m~x,family=binomial('logit'))
## 
## Call:  glm(formula = y/m ~ x, family = binomial("logit"))
## 
## Coefficients:
## (Intercept)            x  
##       -2.50         3.29  
## 
## Degrees of Freedom: 19 Total (i.e. Null);  18 Residual
## Null Deviance:       13.9 
## Residual Deviance: 1.66  AIC: 13.5
data=list(N=n,m=m,x=x,y=y)

mdl=cmdstan_model('./ex6-3-0.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -224.264 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       10      -77.1961    0.00156079   2.52577e-06           1           1       13    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -77.20
##      b0      -2.50
##      b1       3.29
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -78.18 -77.87 1.00 0.70 -80.18 -77.25 1.00      748      767
##      b0    -2.58  -2.56 0.45 0.44  -3.35  -1.88 1.01      526      489
##      b1     3.38   3.37 0.47 0.47   2.62   4.17 1.01      535      552

b0=mcmc$draws('b0') |> as.vector()
b1=mcmc$draws('b1') |> as.vector()

pse=-b0/b1
quantile(pse,probs=c(0.0,0.025,0.05,0.25,0.5,0.75,0.95,0.975,1),na.rm=T)
##    0%  2.5%    5%   25%   50%   75%   95% 97.5%  100% 
## 0.572 0.645 0.665 0.721 0.760 0.802 0.866 0.880 0.964
jnd=(log(0.75/0.25)-b0)/b1-pse
quantile(jnd,probs=c(0.0,0.025,0.05,0.25,0.5,0.75,0.95,0.975,1),na.rm=T)
##    0%  2.5%    5%   25%   50%   75%   95% 97.5%  100% 
## 0.219 0.256 0.263 0.298 0.326 0.359 0.420 0.442 0.547
x1=runif(length(b0),-2,2)
p=1/(1+exp(-(b0+b1*x1)))
pm=1/(1+exp(-(median(b0)+median(b1)*x1)))

qplot(x1,pm,col=I('darkgray'),ylab='p')+
  geom_line(aes(x=x1,p=p),col='red')